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We analyze the synchronization transition for a pair of 
coupled identical Kauffman networks in the chaotic phase. 
The annealed model for Kauffman networks shows that syn- 
chronization appears through a transcritical bifurcation, and 
provides an approximate description for the whole dynam- 
ics of the coupled networks. We show that these analytical 
predictions are in good agreement with numerical results for 
sufficiently large networks, and study finite-size effects in de- 
tail. Preliminary analytical and numerical results for partially 
disordered networks are also presented. 
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I. INTRODUCTION 

Synchronization of coupled elements is a form of col- 
lective evolution present in a variety of complex real sys- 
tems and mathematical models. This class of emergent 
behavior has been observed in biological populations |y , 



chemical reactions 
social phenomena 



neural networks and human 
among other instances. Models 
that account for synchronization consider, for example, 
globally coupled logistic maps Q, chaotic oscillators 
Hamiltonian systems 0j, and formal neural networks jgj. 

In the usual formulation, two identical dynamical sys- 
tems whose individual dynamics is governed by the equa- 
tion w = F(w) are coupled to each other in the form 



Wi,2 = F(wi j2 ) + e (w 2 ,i - wi j2 ) . 



(1) 



where e is the coupling parameter. Full synchronization 
takes place when both systems converge asymptotically 
to a common trajectory, wi(t) = W2(f). When the indi- 
vidual dynamics is chaotic — a particularly relevant case 
in connection with the description of real systems — full 
synchronization occurs above a critical value e c of the 
coupling intensity. This critical point is determined by 
the competition between chaos and coupling, and can be 
calculated in terms of the Lyapunov exponent of the in- 
dividual dynamics ||. 

While globally coupled chaotic elements with a few in- 
ternal variables have been extensively studied, synchro- 
nization of spatially extended systems remains quite un- 
explored. Recently, synchronization has been reported 
for a system consisting of two coupled complex Ginzburg- 
Landau equations Globally coupled neural networks 
H , stochastically coupled cellular automata [ ^0| - [T2|, and 
nonidentical complex Ginzburg-Landau systems ]|13[| , are 



other examples of spatially extended systems that present 
a critical transition to synchronization. 

In this paper we study the synchronization dynamics of 
two coupled identical Kauffman networks, which are dis- 
crete extended dynamical systems with quenched disor- 
der. With respect to previous work on coupled extended 
systems, the interest of Kauffman networks resides in the 
fact that the transition to synchronization admits an an- 
alytical description which results to be in excellent agree- 
ment with numerical simulations. In Sect. [H] we briefly 
review the definition of Kauffman networks and the an- 
nealed model for the calculation of their overlaps. Next, 
in Sect. [II, wc introduce a stochastic coupling mecha- 
nism for Kauffman networks and propose an analytical 
approach in the framework of the annealed model, which 
identifies the transition to synchronization in our system 
as a transcritical bifurcation. Section [V , where we report 



our numerical results, is the core of the present paper. 
There, we study the effects of spurious synchronization 
in finite-size networks, consider the application of noise 
to the system to eliminate such effects, and compare the 
results with the analytical description. Remanent finite- 
size effects are numerically quantified and their analytical 
treatment, which requires a formulation beyond the an- 
nealed model, is outlined. In Sect. [V|we discuss the syn- 
chronization transition in some subclasses of Kauffman 
networks, which may be thought of as interpolations be- 
tween generic Kauffman networks and ordered cellular 
Finally, in Sect. VI, our results are summa- 



automata 
rized and discussed 



II. KAUFFMAN NETWORKS 

Kauffman networks, also known as random Boolean 
networks, were introduced as a model for the problem of 
cell differentiation |lj,|l^]. Since then, they have been 
the object of many studies concerning their properties 

A Kauffman network (KN) is a disordered determinis- 
tic dynamical system. It consists of an iV-site network, 
where each site is connected to K randomly chosen sites. 
The parameter K is known as the connectivity of the net- 
work. We refer to the set of K sites connected to a given 
site as its neighborhood. The state of each site is given by 
a Boolean variable Oi S {0, 1}, and evolves according to 
the inputs coming from the neighbor sites. The evolution 
rule is chosen independently and randomly for each site. 
To each possible configuration of the neighborhood — 
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there are 2 such configurations — an output is assigned, 
namely, 1 with probability p, or with probability 1 — p. 
The parameter p e [0, 1] is known as the bias of the rule. 
Then, for each variable <Ji we have Boolean functions fi 
such that Ui(t + 1) = f t [i^(t)], where v { = (a^,. ..,a iK ) 
is the set of inputs of site i. The state of all sites is 
updated simultaneously according to the corresponding 
functions. We can write an evolution equation for the 
state vector of the network a = (p\, 02, . . . , cr/v), as 

a(t + l)=f[tr(t)], (2) 

with f [tr(t)] = (/a [i/ x (t)] , h M*)] , • • • Jn 

The if connections and the evolution rule of each site 
are chosen at the beginning and kept fixed during the evo- 
lution. Thus, the disorder is quenched and the dynamics 
is deterministic. For a finite number of sites N, the num- 
ber of states in phase space is also finite — it equals 2 N . 
Then, for any initial condition, the system will eventually 
fall into a cycle. 

In the (p, K) parameter space, Kauffman networks 
present phases of frozen and chaotic evolution, separated 
by a critical line. The transition between these phases has 
been extensively studied, and characterized by means of 
several order parameters, such as the Hamming distance 
jl6],[l7| and the stable core size jl8| . In most of this work 
we will deal with the case p = 1/2 and K = 3, which lies 
within the chaotic phase. 

The annealed model (AM) was introduced to study the 
evolution of overlaps between states in KNs 0, 17TI. In 
this model, the K connections . . . , %k\ of each site as 
well as the Boolean functions fi are randomly changed 
at each time step. This means that an entirely differ- 
ent realization of the network is used at each step. Note 
that, while ordinary KNs are deterministic, the annealed 
model works as a probabilistic automaton. The asymp- 
totic periodic behavior of KNs is absent in the annealed 
model. The advantage of this model is that it allows for 
analytical calculations, and it has been shown that its 
predictions are in good agreement with the behavior of 
KNs in the limit JV->oo @. 

Suppose that we have two identical KNs with the same 
connections and rules. We feed them with different ini- 
tial conditions, and let them evolve in time. We define 
the overlap a(t) between the networks as the fraction of 
homologous sites that are in the same state at time t. In 
the AM, it is possible to calculate the time evolution of 
the overlap. At time t+1 the connections and the rules 
fi are reassigned, but the same changes are applied to 
both networks, keeping them identical. The probability 
for a site having all its inputs coming from sites in the 
same state in both networks is a(t) K . At the next time 
step, consequently, such site will be in the same state 
in both networks, no matter the evolution rule chosen 
for it. Thus, there is a fraction a(t) K of homologous sites 
whose state will coincide at t+1. The remaining 1 — a(t) K 
sites still have a probability of overlapping. Even if the 
state of the neighborhoods of a given site are different in 



the two networks, it may happen that the evolution rule 
assigns the same output to them. The probabilities for 
fi iyl) = fi K ? ) = and for /, (uj) = f t (u 2 ) - 1 are, 
respectively, (1 — p) 2 and p 2 . The overlap at time t + 1 
is then 

a(t + 1) = a{t) K + [1 - a(t) K ] [p 2 + (1 - p) 2 ] . (3) 

An alternative way to characterize the difference be- 
tween two networks is the difference automaton, defined 
by 

di(t)=a}(t)®a 2 (t), (4) 

where ® denotes Boolean addition. The density of this 
automaton is given by 

1 N 

i=i 

and coincides with the Hamming distance between the 
networks. Note that D(t) = 1 — a(t) so that, from Eq. 
(§), we get 

D(t+l)=2p(l-p)(l-[l-D(t)] K ). (6) 

The Hamming distance has proven to be a suitable order 
parameter in the study of the synchronization transition 
in coupled elementary cellular automata p0|-|l2|] , where 
the analysis of overlaps between states is a basic tool to 
define the effects of coupling. In the next section, we 
adapt the annealed model to the description of coupled 
KNs. 



III. COUPLED KAUFFMAN NETWORKS 
A. Stochastic coupling 

In order to establish a coupling mechanism between 
two KNs, we first observe that, due to the discrete na- 
ture of KNs, the usual deterministic coupling used for 
maps H cannot be applied here. Consequently, we in- 
troduce a form of stochastic coupling between networks 
as previously done for cellular automata pi , where the 
continuous parameter q that controls the strength of the 
coupling is a probability, as explained in the following. 

The evolution of the coupled system is implemented 
by the successive application of two operators. First, 
the evolution function f is applied to both networks as 
if they were not coupled [see Eq. (|2j)], yielding f (er 1 ) 
and f (cr 2 ). Next, the stochastic coupling operator S is 
applied: 

{a\t + 1), * 2 (t + 1)} = S (f [a l (t)] , f [cr 2 (t)} ) . (7) 

The operator S compares the states of the networks site 
by site. If crj(t) — <r 2 (t) the state of the site is not modi- 
fied. If, on the other hand, <j}{t) ^ <J 2 (t), with probabil- 
ity q the states of the homologous sites in both networks 
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are set to the same value. This value is chosen among 
<j}{t) and of(t), with the same probability 1/2 for each 
instance. With probability 1 — q, even if aj(t) ^ of (t), 
the coupling does not act, leaving the state of that site 
unchanged in both networks. We call q the coupling prob- 
ability. The whole evolution can be formally expressed 



as 




fi [ v i)>f* ("?) 

fi ( u i)ji 

fi Wi ) > fi ("i) 



with probability 1 — q, 
with probability q/2, 
with probability q/2. 

(8) 



We stress that we are dealing with two identical ex- 
tended systems, and that the coupling mechanism con- 
nects homologous elements of these two systems, namely, 
the ith site of network 1 is connected by coupling with the 
ith site of network 2, as schematically illustrated in Fig. 
[l]. The coupling mechanism defined above is symmetric, 
since each network may influence the other. It could also 
be possible to consider a biased, non symmetric coupling, 
in which one network drives the other Jlltl. 



network 1 



coupling 



network 2 



FIG. 1. Schematic representation of the coupling mecha- 
nism. Coupling links homologous elements of two extended 
systems — in this case, two networks. Since coupling is sym- 
metric, each network may act on the other. 

For q = the two networks are uncoupled, and evolve 
independently from each other. For q = 1, in contrast, 
the networks synchronize completely at the first time 
step. From then on, they follow a common trajectory 
in phase space without further intervention of the cou- 
pling mechanism. Our aim in the following is to study 
the behavior of coupled KNs for intermediate values of 
the coupling probability, q e (0, 1). 

In the frozen phase, were no damage spreading takes 
place, an arbitrary small coupling intensity q > leads 
eventually to synchronization. The situation is differ- 
ent in the chaotic phase. There, we find two competing 
driving forces acting on the coupled system, namely, the 
chaotic dynamics which induces the separation between 
two trajectories to grow and the coupling, by which 
the Hamming distance between the networks decreases. 
In this paper, we focus attention on the chaotic phase. 



B. Annealed model for coupled networks 

The annealed model can be used to predict the behav- 
ior of the pair of coupled KNs. We recall that the time 
evolution equation for the Hamming distance in the case 
of two free networks is given by Eq. (^J). Now suppose 
that the Hamming distance of two coupled KNs at time 
t is D(t). The first substep in the dynamics of this sys- 
tem consists of the free evolution of both networks. The 
Hamming distance after this substep, D(t + 6t), is there- 
fore given by Eq. (g). At the second substep, coupling 
acts on the system, and a fraction q of the homologous 
sites that were in different states are assigned the same 
value. This leaves a fraction (1 — q)D(t + 5t) of sites with 
different states in the two networks. Thus, the evolution 
of the Hamming distance for coupled networks is given 
by the map 



D(t + 1) = F [D{t)\ = cp(p, q) 1 - [1 - D(t)} 



l-K 



(9) 



with tp(p, q) = 2p(l - p)(l - q). 

It can be shown that the map 
point D* > for q < q c , with 

q c = l-[2p(l-p)K}- 1 



has a stable fixed 



(10) 



At q = q c the system undergoes a transcritical bifurcation 
p0| , and D* = becomes a stable fixed point for q > q c . 
Thus, within the AM approximation, q c stands for the 
critical coupling at which synchronization sets on. In 
this paper we deal mostly with the case K — 3, for which 
the stable equilibrium D* can be given analytically as a 
function of the coupling probability q. In this case, in 
fact, the map is defined by the cubic function F(x) — 
(p(p, q)x (3 — 3x + x 2 ) . The stable Hamming distance is 



D*(q) 



3 — 1 
2 2 



-3 



1/2 



for q < q c , 



(11) 







for q > q c 



Note that near the critical point, q < q c , this Hamming 
distance is approximately given by 



D*(q)=6p(l-p)\q-q c \. 



(12) 



Therefore, the corresponding critical exponent is equal 
to unity. 

For q ^ q c , the Hamming distance approaches D*(q) 
exponentially in time. For q — q c , on the other hand, 
Eq. (H) can be approximately written, for D(t) — > 0, 
as D(t + 1) = D(t) - (A" - l)D(t) 2 /2. This implies a 
power-law decay for long times, D(t) ~ t -1 . In the next 
section, we compare these analytical results with those 
of extensive numerical simulations. 
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IV. NUMERICAL RESULTS 

We have performed numerical simulations of pairs of 
KNs coupled under the scheme presented above. The 
results reported in this section correspond to the case 
of p — 1/2 and K = 3. We have recorded the time 
evolution of the Hamming distance, performing averages 
over r realizations of iV-site networks. The number of 
realizations is chosen in such a way that, for different 
values of N, rN > 10 6 . In a typical realization, we start 
with two identical networks with different random initial 
conditions. For each realization, new connections and 
local functions are chosen. The networks are allowed to 
evolve freely, without coupling, for a transient time r 
of, typically, 10 3 steps. This is done for the networks 
to reach their asymptotic dynamics before coupling is 
allowed to act. After this transient, we turn the coupling 
on, reset the time to zero, and start measuring D(t). 

In Fig. H we show the time evolution of the Ham- 
ming distance D(t) for different values of the coupling 
parameter q. The values of q have been chosen to dis- 
play the three typical behaviors, namely, synchronization 
for q > q c , critical decay for q w q c , and convergence to 
a finite distance for q < q c . 
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FIG. 2. Hamming distance as a function of time for 10 3 -site 
networks and three coupling probabilities q, averaged over 10 4 
realizations. The dashed line, of slope —1, is to be compared 
with the power-law decay observed near the critical coupling. 

For the present values of p and K, the annealed model 
predicts a critical coupling probability q c — 1/3 [cf. Eq. 
(|l0|)]. It is however clear from Fig. ||that the power-law 
decay in D(t) is observed for a lower coupling, q = 0.29. 
Simulations of the same networks with q = 1/3, on the 
other hand, always lead to synchronization. In fact, the 
annealed model is expected to provide a good approxi- 
mation to our system in the limit N — > oo jlTfl . Figure || 
shows the time evolution of the Hamming distance for a 
fixed coupling probability, q = 0.29, and several values of 
N. The AM prediction from Eq. (|^) is also shown. The 
strong dependence with the size of the network is appar- 



ent. In particular, we find that for this coupling strength 
10 3 -site networks synchronize whereas 10 4 -site networks 
do not. The AM result gives a good description for the 
case of TV = 10 4 . 
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FIG. 3. Hamming distance as a function of time for several 
network sizes N, and fixed coupling probability q — 0.29. 
Averages were done over 10 3 realizations for N < 10 4 , and 
over 10 2 realizations for N = 10 4 . The dashed curve stands 
for the annealed model prediction. 

The dependence of D(t) on the size of the system can 
be partially ascribed to the effect of spurious synchro- 
nization. Because of the discrete nature of KNs, two 
finite-size networks can be brought to the same state by 
a fluctuation caused by the stochastic coupling. In such 
case, the two networks will remain synchronized from 
then on. This event is more frequent for small networks, 
where the relative size of fluctuations increases. Near 
the critical point, furthermore, where the Hamming dis- 
tance vanishes asymptotically, the effect of fluctuations 
is strongly enhanced. The net result of spurious synchro- 
nization is that the effective critical coupling for finite- 
size networks shifts to lower values as N decreases. As a 
consequence, the average Hamming distance in our simu- 
lations may vanish even for coupling probabilities below 
q c , as illustrated in Figs. [| and 

Spurious synchronization can be avoided by adding 
noise to the system. Such strategy has already been 
adopted in this field, specifically, in the study of globally 
coupled chaotic maps |21 22] ], to prevent synchronization 
due to round-off errors in computer simulations. We im- 
plement the addition of noise as a new substep in the 
dynamics of our system. After the evolution and cou- 
pling substeps, we flip the state of each site in one of the 
networks with a small probability rj. Figure ^ illustrates 
the effect of noise in the evolution of D(t) for q = 0.29 
and N = 10 3 . For this coupling intensity, where the con- 
sequences of spurious synchronization are vast, the be- 
havior of the Hamming distance with and without noise 
changes drastically. 
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FIG. 4. Hamming distance as a function of time for 10 -site 
networks with q — 0.29 and two values of the noise intensity 
77, averaged over 10 3 realizations. The effects of spurious syn- 
chronization for rj = are apparent. 

Note that noise eliminates spurious synchronization for 
q < q c , but also prevents the KNs to exactly synchro- 
nize even for q > q c . Therefore, the critical behavior 
that characterizes the synchronization transition in the 
absence of noise disappears as noise is added, and is re- 
covered only for 77 — » (but 77 7^ 0). The effect of noise 
can be straightforwardly incorporated to the AM approx- 
imation. 
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FIG. 5. Analytical results of the annealed model for the 
asymptotic Hamming distance D*, for different noise inten- 
sities rj. The insert shows the deviation D* n — D% from the 
distance in the absence of noise. 

The map that gives the time evolution of the Hamming 
distance is now 

D(t + 1) = (1 - v)F [D(t)} + 77 (1 - F [D(t)}) , (13) 

with F [D(t)] given by Eq. (||). As for the model without 
noise, for K — 3 it is possible to analytically find the 
asymptotic distance D* predicted by Eq. (p~3|) . Figure 



[| shows D* as a function of the coupling probability q 
for various noise intensities. Note the approximation to 
the critical behavior as 77 — > 0. The insert displays the 
difference between D* and the asymptotic distance in the 
absence of noise as a function of q. 
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networks and three values of the coupling probability q, aver- 
aged over 10 2 realizations. The noise intensity is 77 = 10 -4 . 



In Fig. ^, we compare the prediction of Eq. (13) with 
numerical results for the Hamming distance of two cou- 
pled 10 4 -site KNs with noise intensity 7/ = 10 -4 , for three 
values of the coupling intensity. The agreement is excel- 
lent during the transients, but some noticeable discrepan- 
cies persist in the asymptotic value, especially, for q q c . 




1.0 
FIG. 7 



Asymptotic Hamming distance for three network 
sizes N, averaged over 5 x 10 3 time steps and 10 3 realizations. 
The noise intensity is 77 = 10 -4 . The curve corresponds to the 
AM prediction. 



To study such discrepancies in detail, and thus test the 
AM results, we compute from our numerical simulations 
the average asymptotic value of D(t), defined as 
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, t +T 

(D) = -^D(t). 



(14) 



t=t 



Averages are performed over a time span T of 5 x 10 3 
time steps, when the asymptotic regime of the coupled 
system has been reached, i.e. for sufficiently large values 
of to. As above, we choose rj — 10 -4 , and determine 
(£>) for several values of N. Results for N = 2 7 , 2 10 , 
and 2 14 are presented in Fig. 0. The AM prediction for 
this value of r\ is also shown, as a curve. We see that 
the AM results systematically overestimate the values of 
(D) , and that the agreement improves for larger values of 
N. Therefore, even when noise has been added to avoid 
spurious synchronization, finite-size effects persist. 

These remanent finite-size effects are measured by the 
difference SD* — D* — (D) between the AM prediction 



and the numerical average defined in Eq. (|14|). In Fig. |8j 
we plot SD* as a function of N for different coupling in- 
tensities. The insert shows D* and (D) as a function of N 
for the same values of q. There are two well-differentiated 
regimes in the size dependence of SD*. For small N, the 
deviation between the AM estimate and numerical results 
is practically constant. For large values of N the devi- 
ation decreases, seemingly as a power-law, SD* ~ N~ z . 
Least-square fits for N > 10 2 yield z = 1.1 ± 0.2 for the 
exponent. 



1 N 



(15) 



i=l 



[cf. Eq. (|J)], whereas the distance between two states er 
and er' reads 



1 



A' 



i=i 



a 'i) 2 



(16) 



For fixed p, K and N, a realization R of the network is 
defined by the connections and the local rules. We define 
Q R as the set of all the states visited by the KN for this 
realization, at asymptotically large times and from all 
the possible initial conditions. In other words, the set 
£} R contains all the states that belong to the limit cycles 
of the dynamics. It is possible to introduce a probability 
distribution Vr (<t) over £l R , given by the frequency with 
which a given state er is visited at asymptotically large 
times averaged over all initial conditions. Averages over 
fl R will be computed with this distribution. For instance, 
the average of the density p is 



d R = (p(tr))a R = v n( cr )p( cr )- 

The average of the distance between two states, Eq. 
can be written as 



(17) 
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FIG. 8. Difference between numerical results and the AM 
prediction for the asymptotic Hamming distance as a function 
of size, for four coupling probabilities q. The dashed lines have 
slope —1. The insert shows the asymptotic Hamming distance 
in semilogarithmic scale, to appreciate the behavior for small 
values of N. 

To find an explanation for these finite size effects, it 
is necessary to go beyond the annealed approximation. 
In the following, we outline an approach to the calcu- 
lation of the Hamming distance based on statistical av- 
erages over states of the KNs. The density of a state 
er = (a'i, . . . , ctjv) of a single KN is defined as 



(D(a,a'))a R = 2d R (l 



9 N 

d *>-|£e 9 

i=l 



l > 



(18) 



where & = ((Ti)n R - d R = (a' i )a R - d R . Here, we have 
assumed that the occurrence of states er and er' are prob- 
abilistically uncor related. 

The quantities measure, for each realization of the 
network, the average deviation of the state of each site 
from the average density d R over the whole network. Let 
us introduce the distribution r#(£) as the fraction of sites 
in the KN with deviation £ for a specific realization R. 
Unfortunately, the explicit form of Fr(£) is not known. 
It is however known that this is a nontrivial distribution, 
in particular, due to existence of the so-called stable core 
The stable core is a set of sites that have always 
the same asymptotic, fixed states, irrespectively of the 



initial condition. For these sites, (<Ji 



or 1, so that 



the deviations £j adopt their extremal values, = — d R 
or 1 — d R) respectively. Using the distribution T R (£) to 
replace the sum in Eq. ([Of) by an integral, we have 



(D(<r, tr')) aH = 2d R (l - d R ) - 2 / r R (£)fd£ 



i-df 



(19) 

Now, analytical results on the size of the stable core |l§| ] 
suggest that, as N — > oo, Tr(£) approaches an asymp- 
totic profile r(£) which depends onp and K, but becomes 
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independent of the specific realization of the network. 
For large sizes, we may assume a dependence of the form 

r*(0 « r(0 - a^ 1 ^), (20) 

where T' R is the first (analytical) correction due to finite 
sizes. Within these assumptions, Eq. ( fL9| ) takes the form 

(D(a, <r'))n R = D - | £ 1^(0^^, (21) 
where 

D = 2d fl (l - d«) - 2 / * r(0f 2 dC (22) 

is the asymptotic distance for A?~ — > oo. 

We now associate (D(cr, <t'))q r with the distance be- 
tween the states of two coupled KNs at a given (long) 
time. Indeed, in our system both networks have the same 
connections and rules, and correspond therefore to the 
same realization R of the network. According to the def- 
inition ([IJ) and (|5|), after an average over realizations of 
the network for fixed p, K, and N is performed, the dis- 
tance (D((T, it'))q„ coincides with the Hamming distance 
(D), Eq. (p|), considered above. In Eq. (|l|), thus, D 
should correspond to the Hamming distance predicted 
by the AM approximation, valid for N oo, and the 
first correction to the AM estimate is given by the ad- 
ditional term. As found from our simulations, this term 
depends on the network size as N' 1 . The ansatz (EG) 
is therefore supported by numerical results. Note that 
these conclusions are independent of the strength of cou- 
pling, measured by the probability q, since the only effect 
of the coupling dynamics in connection with the above 
analysis is to change the set of asymptotic states fi^ and 
the profile of the distribution I\r(£). 

V. DISORDERED CELLULAR AUTOMATA AND 
PARTIALLY ORDERED NETWORKS 

According to the results reported in the previous sec- 
tions, the synchronization transition of stochastically 
coupled KNs on the one hand, and of stochastically cou- 
pled elementary cellular automata (ECA, p^|) on the 
other, are qualitatively different. Namely, they belong 
to different universality classes. Whereas we have found 
that synchronization in KNs appears through a transcrit- 
ical bifurcation — with a critical exponent equal to unity 
for the asymptotic Hamming distance — the correspond- 
ing transition in ECA has been shown to exhibit nontriv- 
ial exponents |HJ , which seem to suggest that it belongs 
to the universality class of directed percolation jl2). It 
is therefore relevant to study a third class of systems, 
intermediate between ECA and generic KNs. 

In a generic KN there are two sources of disorder. We 
have the network topology, determined by the random 



choice of connections, and the local evolution rules, which 
are also chosen at random. On the other hand, in cellu- 
lar automata both the topology and the dynamical rules 
are fully homogeneous. ECA can indeed be interpreted 
as a very special subclass of KNs with K = 3, where 
the choice of dynamical rules and connections is deter- 
ministic. In order to distinguish between the effects of 
disorder in the topology and in the dynamics on the syn- 
chronization transition, we consider now the subclass of 
KNs where the connections are still chosen at random 
but the local rule is the same for all sites. We refer to 
these networks as disordered cellular automata (DC A). 

We focus here on the elementary rule 22, which exhibits 
chaotic evolution . This rule is defined by the Boolean 
function /({0,0,1}) = /({(), 1,0}) - /({1,0,0}) = 1 and 
/ = for the remaining five possible neighborhoods (see 
Table |). Thus, the bias for rule 22 is p = 3/8. This DCA, 
however, cannot be thought of as a KN with p — 3/8. In 
a generic KN with this bias, in fact, the local evolution 
could be substantially different from the behavior of rule 
22. In particular, the dynamics at some sites may be 
governed by nonchaotic rules, giving rise to sensible dif- 
ferences in the global behavior. This is clearly illustrated, 
for instance, by a measurement of the asymptotic den- 
sity of a rule-22 DCA, which yields d ~ 0.423, instead of 
d = p = 0.375. 

The formulation of an annealed model of DCA requires 
a different approach, in order to account for the homo- 
geneity of the dynamical rules. We define the AM by 
reassigning all the connections at each time step, but 
keeping the functions fi = f fixed. Suppose that we 
have two networks with overlap a(t). The probability 
for a site to have exactly K — I of its reassigned inputs 
coming from homologous sites in the same state is 

Pi{t)=(^ja{t) K - l [l-a{t)} 1 . (23) 

We introduce the quantity Ai as the probability for two 
homologous sites having all but I inputs coming from ho- 
mologous sites in the same state to give the same output. 
The overlap at time t + 1 is then given by 

K 

a(t + l) = ^2Aip,(t). (24) 

1=0 

The quantities Ai depend on the evolution rule, and 
their value is fixed. They can be evaluated within some 
approximations, as shown in the Appendix. Note that 
Aq = 1 because, no matter the rule, if the inputs are all 
equal the outputs will coincide. In the annealed model 
for KNs we had A = 1 and Ai — p 2 + (1 — p) 2 for 
I = 1,...,K. For K — 3 the map for the Hamming 
distance can be casted into the form 

D(t + 1) = BiD(t) + B 2 D(t) 2 + B 3 D{tf, (25) 

with B x = 3(1-Ai), B 2 = 3 (2A X - A 2 - 1), B 3 = 
—3Ai + 3A 2 — As + 1. Coupling enters then the for- 
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mulation exactly as in Eq. (g), as an additional factor 
1 — q in the evolution of D(t). 

In Fig. ^, the curve stands for the asymptotic value 
D* as a function of the coupling probability predicted 
from Eq. (|25|) for rule 22. The prediction is qualitatively 
similar to that for KNs, in particular, in the region close 
to the synchronization transition. Numerical results on 
DCA with rule 22 for N = 2 11 are also shown in Fig. 
|9|. To avoid spurious synchronization, a small amount 
of noise, r\ = 10~ 4 , has been added. The agreement 
with the AM is reasonably good, though some systematic 
deviations are clearly visible in the zone of the transition. 
As before, these deviations may be attributed to finite- 
size effects. 




FIG. 9. Numerical results for file asymptotic Hamming dis- 
tance in coupled 2 -site disordered cellular automata with 
the evolution rule 22, for different amounts of disorder. The 
curve stands for the annealed model prediction for DCA. 



Finally, as an interpolation between DCA and ECA 
we consider partially ordered topologies (POCA), con- 
structed in the following way. We start with an iV-site 
ECA, which consists of a one-dimensional array where 
each site is connected to itself and to its two nearest 
neighbors. Then, for each site, each neighbor is redrawn 
at random from the whole network with probability 9. 
With the complementary probability, 1 — 6*, the neigh- 
bor is preserved. The resulting topological structure of 
connections is analogous to that of small-world networks 
p4j . For 9 = 1, we recover the DCA networks discussed 
above. Results of numerical simulations on POCA with 
rule 22 are shown in Fig. || for 9 = 0.1 and 9 = 0.3. 
Data for ECA with rule 22 (9 = 0) are also shown. We 
see that, in spite of the relatively low values of 9, the 
asymptotic Hamming distance of POCA depends on the 
coupling probability q in a way qualitatively similar to 
that of KNs. In particular, it does not exhibit the abrupt 
dependence on q observed in ECA near the critical point. 
This would be in agreement with the crossover scenario 
in small- world networks |25j , where it is known that even 
small amounts of disorder induce behaviors which already 



resemble that of fully random structures. 

VI. SUMMARY AND DISCUSSION 

We have studied the behavior of two Kauffman net- 
works interacting through a form of symmetric stochastic 
coupling. As for many other localized or extended, ran- 
dom or deterministic, coupled dynamical systems P-p^[, 
we have found that coupled Kauffman networks can syn- 
chronize their evolution if coupling is strong enough. In 
our case, there is a critical value of the coupling proba- 
bility q beyond which the two networks converge to the 
same trajectory as time elapses. 

In contrast with the situation encountered for other ex- 
tended systems |q-|l3|], however, for Kauffman networks 
it has been possible to give an analytical description of 
the synchronization transition, in excellent agreement 
with numerical results for large-size systems. This for- 
mulation is provided by an extension of the so-called an- 
nealed model to the system of coupled networks. 
The model gives the evolution of the overlap between 
the two networks — or, in other words, of their Hamming 
distance — and makes it possible to evaluate its asymp- 
totic value. The asymptotic Hamming distance D* is 
used as an order parameter for the synchronization tran- 
sition. For coupled Kauffman networks in the chaotic 
phase, the annealed model predicts the existence of a 
critical coupling probability q Cl such that D* is finite for 
q < q c and vanishes for q > q c . At the critical point, 
D* ~ | (7 — <7 C | , and the transition has the character of a 
transcritical bifurcation. 

The transition predicted by the annealed model is qual- 
itatively different from that observed in the deterministic 
version of this kind of systems, more specifically, in cellu- 
lar automata. For elementary cellular automata, in fact, 
the critical behavior of the Hamming distance exhibits 
nontrivial exponents pfj| . It has been suggested that, at 
least for some evolution rules, the synchronization transi- 
tion in cellular automata belongs to the class of directed 
percolation Jl^] . On the other hand, by analogy with the 
problem of damage spreading, synchronization in Kauff- 
man networks could belong to the class of directed per- 
colation in disordered systems ptsfl . 

Results from extensive numerical simulations of rela- 
tively large networks (10 4 sites), performed for Kauff- 
man networks in the chaotic regime (p = 1/2, K = 3), 
are in good agreement with the predictions of the an- 
nealed model both in the temporal evolution and in the 
asymptotic behavior of the system. However, for small 
networks some systematic departures from the analytical 
results are apparent. These deviations can be partially 
explained taking into account the occasional events of 
spurious synchronization for q < q c , due to the effect of 
suitably large fluctuations on our discrete finite-size sys- 
tem. The effect of these fluctuations becomes more im- 
portant as the coupling probability approaches its critical 
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value. In any case, spurious synchronization can be suc- 
cessfully eliminated by adding a small amount of noise to 
the dynamics. Moreover, noise can also be encompassed 
into the annealed model. We have introduced noise in the 
simulations and observed that remanent finite-size effects 
persist. These must now be ascribed to the failure of the 
annealed model in describing finite networks. By means 
of a more detailed statistical description of the overlaps 
between networks, in fact, we have been able to account 
for such remanent effects, also explaining the dependence 
of the deviation from the annealed model on the network 
size. 

Finally, we have presented preliminary numerical re- 
sults on disordered and partially ordered stochastically 
coupled cellular automata. These systems can be seen 
as providing a kind of interpolation between Kauffman 
networks, with their completely random connections and 
dynamical rules, and cellular automata, which are fully 
ordered. In the case of disordered networks with the same 
evolution rule on all the sites, it is possible to extend the 
annealed model, which predicts the same class of synchro- 
nization transition as for Kauffman networks. Numerical 
simulations agree with these predictions. Increasing the 
order in the connections by means of a scheme analo- 
gous to the construction of a small- world network 
we have also considered partially disordered structures. 
Even for small amounts of structural disorder, the behav- 
ior associated with the synchronization transition seems 
to resemble that of Kauffman networks more than that 
of cellular automata. This leads us to conjecture that 
the synchronization transition of partially disordered au- 
tomata is in the same universality class as for coupled 
Kauffman networks. Nevertheless, further extensive sim- 
ulations and a careful determination of the critical ex- 
ponents are necessary to draw any conclusions on this 
point. 
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APPENDIX: COEFFICIENTS FOR THE 
ANNEALED MODEL 



we neglect the correlations between neighborhood fre- 
quencies, associated with the spatial patterns generated 
by the dynamics. For K = 3, there are eight possible 
neighborhoods, which we label from to 7 as shown in 
Table |. Within the above assumption, the frequency pi 
of each neighborhood % can be estimated in terms of the 
density d as 



Pi = P2 = Pa = d(l - d) 2 , 
P3 =P5=Pe = d 2 (l - d), 
P7 = d 3 - 



(Al) 



Moreover, we note that the density d can in turn be writ- 
ten in terms of the frequencies pi of neighborhoods with 
nonzero output. For rule 22 these neighborhoods are 
{1,0,0}, {0,1,0} and {0,0,1} (see Table |), so that we 
have 

d = pi+p 2 +Pi- (A2) 

Combining Eqs. Q and (|A|) yields d = 1 - \/3/3 « 
0.423, which agrees with the numerical measurement re- 
ported in Sect. The corresponding values of the fre- 
quencies pi are shown in Table [|. They are in very good 
agreement with the values obtained from numerical simu- 
lations, also shown in the table. This supports our above 
assumption of uncorrelated neighborhood frequencies. 

Once the frequencies pi are known, we are able to cal- 
culate the coefficients Ai. As a specific example, we dis- 
cuss A3. By definition, this is the probability for two 
homologous sites having all but 3 inputs coming from 
homologous sites in the same state to give the same out- 
put. Thus, we are in the case where the two neighbor- 
hoods have all the homologous sites in different states. 
The pairs of neighborhoods that satisfy this condition 
are {0,7}, {1,6}, {2,5}, and {3,4}. Among them, how- 
ever, only the pair {0,7} has the same output for both 
neighborhoods (see Table |) . The coefficient A3 is there- 
fore given by 



P0P7 



P0P7 + P1P6 + P2P5 + P3P4, 



(A3) 



The computation of the other coefficients is accomplished 
in a similar way, yielding A x = 4(19 - 8\/3)/169 and 
A 2 = (9 + V3)/13. Moreover, A = 1. 



In this Appendix we illustrate the computation of the 
coefficients Ai in Eq. (|24|) with an explicit example. We 
recall that Ai is defined as the probability for two homol- 
ogous sites having all but / inputs coming from homolo- 
gous sites in the same state to give the same output. 

Let us assume that the frequency with which a given 
neighborhood appears in a state of the whole network is 
fully determined by the density d of the state. Namely, 
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label 


neighborhood 


output 


^numerical 


/^analytical 





{0,0,0} 





0.19350 


0.19294 


1 


{0,0,1} 


1 


0.14075 


0.14096 


2 


{0,1,0} 


1 


0.14071 


0.14096 


3 


{0,1,1} 





0.10287 


0.10298 


4 


{1,0,0} 


1 


0.14079 


0.14096 


5 


{1,0,1} 





0.10290 


0.10298 


6 


{1,1,0} 





0.10289 


0.10298 


7 


{1,1,1} 





0.07559 


0.07524 



TABLE I. Cellular automata neighborhoods and their 
outputs for the evolution rule 22. The frequencies p for 
each neighborhood obtained from numerical results and from 
the analytical approximation used in the Appendix are also 
quoted. 
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